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This paper addresses the problem of recovering the spatial mar- 
ket potential of a retail product from spatially distributed sales data. 
In order to tackle the problem in a general way, the concept of spatial 
potential is introduced. The potential is concurrently measured at dif- 
ferent spatial locations and the measurements are analyzed in order 
to recover the spatial potential. The measuring instruments used to 
collect the data interact with each other, that is, the measurement at 
a given spatial location is affected by the concurrent measurements 
at other locations. An approach based on a novel geostatistical model 
is developed. In particular, the model is able to handle both the mea- 
suring instrument interaction and the missing data. A model estima- 
tion procedure based on the expectation-maximization algorithm is 
provided as well as standard inferential tools. The model is applied 
to the estimation of the spatial market potential of a newspaper for 
the city of Bergamo, Italy. The estimated spatial market potential is 
eventually analyzed in order to identify the areas with the highest po- 
tential, to identify the areas where it is profitable to open additional 
newsstands and to evaluate the newspaper total market volume of 
the city. 



1. Introduction. The market potential of a given retail product is the 
expected sales volume when the product is marketed. The spatial market 
potential is the spatial distribution of the market potential over a trading 
area. Sales are expected to be high if a store is opened at a spatial location 
characterized by a high spatial market potential while they are expected to 
be low if the spatial location has a low spatial market potential. 

With the goal to increase and to maximize the sales volume, a key issue 
is how to evaluate the spatial market potential. In this paper, it is assumed 
that the product is already marketed and that the sales data of spatially 
distributed stores are available. Thus, the aim is to estimate the spatial 
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market potential by considering the sales data, the spatial characteristics of 
the trading area and the spatial interaction between the stores. The stores 
interact in the sense that the sales volume of each store is affected by the 
presence of all the others. As a consequence, the spatial market potential 
cannot be estimated ignoring the interaction. 

For all purposes and intents, the spatial market potential can be regarded 
as a spatial surface as it is well defined for all the spatial locations of the 
trading area. Taking a statistical perspective, the spatial market potential is 
considered as a spatially continuous random field and the estimation of the 
spatial market potential is obtained through the estimation of the realization 
of the random field. Although well understood, however, no attempt has ever 
been made to address the problem following a geostatistical approach. 

The estimation of a spatial market potential is an instance of the more 
general problem of recovering the realization of a spatially continuous ran- 
dom field in the case of interacting measuring instruments. The instruments 
interact in the sense that the measurement at a given spatial location is 
affected by the concurrent measurements at nearby locations. 

A novel model able to handle both the interaction between the measuring 
instruments and the missing data is proposed. A case study is presented, in 
which the sales data of spatially distributed newsstands are used to estimate 
the spatial market potential of an economic daily newspaper for the city of 
Bergamo, Italy. The aim of the study is threefold: to identify the areas with 
the highest market potential, to identify the areas where it is profitable to 
open additional newsstands and to estimate the total market volume of the 
city with respect to the newspaper considered. 

The rest of the paper is organized as follows: Section 2 provides the back- 
ground and motivation for this work. Section 3 introduces a novel geostatis- 
tical model for the analysis of spatial point data in the case of interaction 
between the measuring instruments. Model estimation and inference are 
discussed in Section 4. Section 5 presents the case study while Section 6 
provides conclusions. The technical aspects related to the model estimation 
are reported in the Appendices. 

2. Background. In this section, the spatial market potential estima- 
tion problem is discussed in terms of both the current state of the art and 
the available statistical methods. It turns out that the estimation of a spa- 
tial market potential from sales data received little attention in the past 
and that the classic geostatistical approach cannot be adopted to solve the 
problem. 
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2.1. Spatial market potential estimation. The problem of estimating the 
market potential of a retail product is not new in the literature. The state 
of the art is represented by the so-called spatial interaction models which 
describe the market potential in terms of flows between a set of origins (the 
customers) and a set of destinations (the stores). The interested reader is 
referred to the seminal papers of Reilly and Huff (Reilly (1931), Huff (1964)) 
and to the more recent literature (see, for instance, Davis (2006), Cliquet 
(2006), de Grange, Ibeas and Gonzalez (2009) and Fischer and Wang (2011)) 
The spatial interaction models focus on the utility that consumers obtain 
from buying a retail product at a specific store. The utility is often a func- 
tion of the attributes of the product, the attributes of the store and some 
attributes of separation such as the geographic distance, the transport cost 
and the transport time. 

The main drawback of the spatial interaction models is that the spatial 
market potential is not explicitly modeled as a regionalized variable and 
it is defined only at the spatial location of the stores. This is in contrast 
with the concept of spatial market potential adopted in this paper, which 
is supposed to exists beyond the existence of the stores. Indeed, the market 
potential of a product at a given location in space can also be considered as 
the willingness of the consumers to reach that specific location in order to 
buy the product. The attributes of the store (including the price at which 
the store sells the product) may affect the way the spatial market potential 
is observed but the spatial market potential is not driven by the stores. 

The spatial interaction models literature also lacks of methods for esti- 
mating the model output uncertainty. This represents a critical issue as, in 
practical applications, the available data are usually limited in number and 
the reliability of the model output must be provided. 

Modelling the market potential as regionalized variable, and in particular 
as a spatially continuous random field, allows to answer new and interesting 
questions. Denoting q{s) the spatial market potential at the generic spatial 
location s, the company that owns the stores may be interested in estimating 
g(s) for each point of the trading area V. For instance, the company may 
want to locate the maxima of q{s) to be sure that it has a store near that 
location. If, instead, the company wants to open a new store, then it may 
want to evaluate the market potential conditioned on the presence of the 
actual stores. Moreover, the company may want to pursue both of the goals 
even if the sales data of some stores are missing. In this sense, when a latent 
spatial market potential has to be assessed and the uncertainty information 
must be provided, the geostatistical approach seems to be more appropriate 
than any approach based on the spatial interaction models. 
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2.2. Geo statistical modeling. The problem of estimating a spatially con- 
tinuous random field from measurements collected at a finite number of 
locations in space is usually solved by considering geostatistical models and 
kriging techniques (see Cressie and Wikle (2011)). 

The simplest geostatistical model described in Diggle and Ribeiro (2007), 
for instance, assumes an underlying stationary Gaussian random field w{s) 
and observation y(sj) which are realizations of conditionally mutually in- 
dependent random variables Y (sj) conditionally normally distributed with 
mean E (Y (sj) | w (■)) = w (sj) and variance cr^. In the spatial market poten- 
tial case, however, due to the interaction between the stores, the conditional 
mutual independence of the random variables Y (sj) is not met and, in gen- 
eral, E (Y (si) I (?(•)) 7^ ^(si). A geostatistical approach for spatial interac- 
tion data can be found in Banerjee, Gelfand and Polasek (2000) though the 
interaction is defined in terms of flows between destinations and origins (cfr. 
the previous paragraph) and the approach is not suitable for the problem 
addressed in this paper, where a "diffuse" origin is considered. 

3. The geostatistical potential model. 

3.1. Introduction. Before developing a suitable geostatistical model, the 
problem at hand is restated and generalized in the following way. 

Let q{s) be a spatial random field defined over the region of space D cM?. 
The random field q{s) is called here potential, though the term does not 
refer to any particular property of the field. The random field is concurrently 
measured at the set of spatial locations S = {si, sat} and the observations 
y (5) are collected (possibly with missing data). The concurrency of the 
measurements is a key aspect in the sense that, in general, the observations 
y (S) collected in the case of non-concurrent measurements differ from y (S). 

The observations y (S) are supposed to be realizations of random variables 
Y (sj) conditionally normally distributed with conditional mean 

E{Y{si) \q{-),S) = h{qis,);S); Si e S,i = 1, N 

and conditional variance cr^, where /i is a function modeling the interaction 
between the measuring instruments. 

The interaction between the measuring instruments is said to be an ab- 
sorption interaction if, for each s G D, /i (g (sj) ; 5) < q{s). The interaction 
is such that, for each Si £ V, h (q (si) ; 5) = g (si) iff 5 = {si}, that is, the 
conditional mean of Y (si) is equal to q (si) if si is the only spatial location 
of T> where q is measured. Ultimately, it can be stated that the act of mea- 
suring the potential q at a given s alters the (concurrent) measurements at 
other spatial locations. 
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At this point, the fohowing distinction can be made: the potential q{s) is 
the expected observed value when q is measured only at the spatial location 
s G P. On the other hand, the conditional potential q{s; S) is the expected 
observed value when q is measured at the spatial location s G D given that 
it is concurrently measured at the set of locations S = {si, ...,S7v}, Sj G V, 



In the next paragraph, the way the potential and the conditional potential 
are modeled and estimated is discussed. 

3.2. Model definition. The geostatistical potential model (GPM) is intro- 
duced here as the main statistical tool for the analysis of spatial data arising 
from concurrent measurements in the presence of interaction between the 
measuring instruments. In its general form, the GPM is described by the 
following hierarchy of equations: 



At the first stage of (3.1), y{s;S) is the measured conditional potential at 
the spatial location s while : M x D x E — y M is the interaction function 
which is parametrized by the parameter vector i9. The set § is the set of all 
finite spatial point patterns over V including the non-simple patterns (i.e. 
patterns with overlapping points). At the second stage, e(s) represents an 
error component which is assumed to be i.i.d. N{0,a'^) and is supposed to 
capture both the measuring error and the model error. Finally, at the third 
stage, the potential q{s) is modeled by three summands, where /i is the mean, 
x(s) is a vector of covariates, /3 is the vector of coefficient, w^s) is a zero-mean 
latent Gaussian process and 7 is a scale parameter. The covariance function 
of w{s) is cov {w{s),'w{s')) = p0{s,s'), with p0(s, s') a valid correlation 
function parametrized by the vector 6. The model parameter vector is ^ = 
/3', C7^, 7, 0', 1?') and it completely characterizes the GPM. 

Note that, for the reasons discussed later on in the paper, it is assumed 
that, conditionally to the observed covariates x(s), the observed y{s;S) is 
not preferentially sampled (see Diggle, Menezes and Su (2010)) with respect 
to the latent variable w. As a consequence, the set of spatial locations S is 
treated as a constant rather than as the realization of a spatial point process 
the local density of which is driven by w. 

In order to have a better insight into the role of the interaction function 
/i^, the following family of interaction functions is adopted: 



iV > 1. 



(3.1) 



y(s;5) 

n(s) 
q{s) 



(n(s) ,s,5) 
g(s) +e(s) 
ji + x(s)/3 + 7'u;(s) 
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(3.2) /i^ (n(s) ,s,5) = ^(s) • I 1 + ^ (s,s') j = u (s) • ^(^(s; 5) 

where (s,s') : x — > M"*" is a generic non-negative binary function. 

The function (s, s') can be any continuous function but, for practical 
appUcations, it should be monotonically decreasing with respect to distance. 
For instance, 

(3.3) (s,s') = (||s - s'll) = exp 

where ||-|| is the Euclidean distance and tD = {(f), a)' is the function parameter 
vector. In equation (3.3), (j) defines the strength of the interaction while a > 
is a shape parameter. 
Note that 

y(s;<S) = u (s) ■ g^{s; S) 

(3.4) = q{s) ■ g^{s;S) + e{s) ■ g^{s;S) 

= q{s;S) + e{s;S) 

namely the observed potential is equal to the conditional potential q{s; S) 
plus a transformation of the error e(s). In particular, the second line of 
equation (3.4) follows directly from the second stage of model (3.1) while 
the third line is the second line rewritten in a more compact notation. 

The term g^(s;S) is the key element of the interaction function and it 
deserves more explanation. If, as an example, the function (3.3) is considered 
and 5 = 0, namely if there are no measuring instruments, then (/^(s; 5) = 1 
since the summand in equation (3.2) cannot be evaluated and it is equal to 
zero by definition. When a measuring instrument is added, S = {si}, the 
potential at s is a function of the distance between s and si. In particular, 
if s = Si then g^{s;S) = 0.5. On the contrary, when ||s — Si|| — )• oo then 
5^(s; 5) ^> 1. This reflects the fact that the action of absorbing the potential 
at site si influences the measure at the site s. It is worth noting that s and 
si are exchangeable in the sense that absorbing and measuring the potential 
are equivalent actions and that the potential cannot be measured without 
being absorbed. 

In this work, the measuring instruments are supposed to be equally- 
effective, that is, g{si;{sj}) = g{sj;{si}) for all ||sj — Sj||. The property 
of equally-effectiveness is satisfied if the binary function (s, s') is com- 
mutative^, which is the case of the function (3.3). In practical applications, 

^The binary function / is commutative if f{x,y) = f{y,x). 
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the property may not be satisfied in the sense that a measuring instrument 
might be more effective in absorbing the potential than a second instrument 
close to it. Suppose equally-effective measuring instruments, however, sim- 
plify the model and any discrepancy from it is accounted by the error term e. 
Note that the measure of effectiveness is strictly related to the the measure 
of attractiveness of the spatial behavior of consumers models typical of the 
geomarketing literature (see Cliquet (2006)). 

To better understand the notions of potential and conditional potential, 
the following example is considered. Suppose that four measuring instru- 
ments are located at si = (0.2,0.2), S2 = (0.2,0.8), S3 = (0.8,0.2) and 
S4 = (0.8,0.8), Si e P = [0, 1] X [0, 1], i = 1, ...,4. The GPM considered is 

yis;S) = u{s) ■ g^{s;S) = q{s-S) 

(3.5) n(s) = q{s) 
q{s) = w{s) 

namely it is supposed that the conditional potential q{s; S) is observed with- 
out error. Furthermore, suppose that 

(3.6) P0{s,s') = p0{\\s 

(3.7) /^(s,s') = /^(||s 

and that y(sj; 5 \ Sj) = 10.^ The estimated potential and conditional poten- 
tial are reported in the left and in the right part of Figure 1, respectively. 
Regarding the potential, its value at the measuring instrument locations 
is equal to 13.2 > 10. Each measuring instrument measures a potential 
equal to 10 since a fraction of it is absorbed by the remaining measuring 
instruments. Indeed, the potential (?(sj) = 13.2 would be measured by the 
single measuring instrument if the other instruments were not present. The 
conditional potential, as expected, has its lowest value at the measuring in- 
strument locations and represents the potential that would be observed by 
a fifth measuring instrument if placed at the generic s. 

4. Parameter estimation and inference. Let y = y (5) be the A^x 1 

vector of data collected at the sampling sites S. The measurement equation 
for the vector y is 

(4.1) y = G(l^ + X/3 + 7w + e) 

■^Note that, in general, y{s;S) should be simulated following equation (4.12). In this 
case, in order to better appreciate the role of the interaction function, y{si;S \ Si) is 
supposed to be equal for all the locations. 
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Fig 1. (left) potential q{s) ; (right) conditional potential q{s; S) . 



where 1 is the x 1 vector of ones, X = X (S) is the N x b matrix of 
covariates, w = w (S) is the latent Gaussian process at S with variance- 
covariance matrix Sw = (5, 6) and e = e (S) is the measurement error 



at S with diagonal variance-covariance matrix Sl^ = a^I^. Finally, G = 
(S) is the N x N diagonal matrix whose diagonal vector is 

g= {gi}{si;S\si),...,gi}{sN;S\sN)) 

Furthermore, suppose that S is partitioned as 5*^^-*}, where S^^^ 

is the set of sites where the data are available and 5^^^ is the set of sites 
where the data are missing. According to this, the vector y is partitioned as 
y* = (y*-^-*, y^^^) , where y*-^-* = Ly is the sub- vector of the non- missing data 
and L is the appropriate elimination matrix. The vector y* is a permutation 
of y and y = Dy*, with D the proper commutation matrix. The partitioned 
measurement equation becomes 

yW = g(^) ( l(*V + X(*)/3 + 7wW + eW) ; i = i, 2 



and the variance-covariance matrix of the permuted errors is conformably 
partitioned as 

Rii Ri2 



Var 



R,21 R'22 



In the sequel, given b a generic vector and B a generic matrix, b^^^ and B^^^ 
will stand for Lb and LBL', respectively, bearing in mind that, in general 
LB-^L' / (LBL')"V 

Given the data vector y and considering the GPM, the following inferen- 
tial problems are of interest: 
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1. to provide an estimate of the model parameter vector 

2. to provide confidence intervals for the elements of ^; 

3. to estimate the potential q (s) over the region V and its uncertainty; 

4. to estimate the conditional potential q{s; S) over the region V and its 
uncertainty; 

5. to evaluate the expected total potential measured by a maximum of 
measuring instruments. 

4.1. Parameter estimation. Problem 1. is tackled here following the max- 
imum likelihood (ML) approach. Being w (s) a latent process and due to pos- 
sible missing data, the expectation-maximization (EM) algorithm is adopted 
to find the ML estimate I' of 

The EM algorithm is based on the complete-data likelihood function 
Lij, (y, w) and it provides an iterative procedure to update the model param- 
eter estimate from ^(*^) to ^(^^+1) until convergence (see McLachlan and Krishnan 
(2008)). In particular, for each iteration of the algorithm, the E-step com- 
putes the conditional expectation 



Qh,¥>'A =E.,,, (y,w 



.(1) 



while, at the M-step, the following maximization is carried out: 

^ argmaxQ (^^, i/'^^^ 

which is equivalent to solve the equation 



(4.2) 



5^' 



0. 



Considering the approach described in Fasso and Finazzi (2011), the fol- 
lowing closed form updating formulas have been derived: 



(4.3) 
(4.4) 

(4.5) 
(4.6) 



tr 



(eW + ^^Wld)) (1(1))'' 



N -Nrr 



1 / e«-(e«)' + (7W)'A« 



N 



tr 



R 



■22 



(eW+^WwW) (wW)' 



tr w(i) (w(i))' + A(i) 
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where e(i) = (G^^)) ^ y(i)-//('=)l(i) -X(i)/3('=) -7('=)w(i), iV^ is the num- 
ber of missing data in y and 

(4.7) w = (w I y(i)) 

(4.8) A = Var^^u) (w | y(i)) 

are the estimated latent variable and the estimation variance, respectively. 
The evaluation of (4.7) and (4.8) is reported in Appendix A. 

The remaining model parameters can be updated by numerical optimiza- 
tion solving ( ^(^^+1), ijC^+i) J = argmaxQ ( ^'^^^ ) . If both the correlation 

function pg and the interaction function have analytical form of the first 
and second derivative with respect to 9 and i9, respectively, both the pa- 
rameters can be updated by adapting the algorithm given in Xu and Wikle 
(2007). 

Before concluding the paragraph, a point that is worth mentioning is how 
the preferential sampling problem can affect the model parameter estimation 
in the case of the GPM. The spatial data y (5) are preferentially sampled 
with respect to the potential q{s) if the spatial pattern of S is not indepen- 
dent of g(s). In practice, the spatial density of S can be higher at the spatial 
locations where g(s) is known or expected to be high. 

As discussed in Diggle, Menezes and Su (2010), if the data are preferen- 
tially sampled and the issue is not addressed, then the estimation of the 
parameter 9 related to the latent variable w{s) is generally biased. With re- 
spect to the GPM, however, the model considered in Diggle, Menezes and Su 
(2010) does not include covariates. If the data are preferentially sampled with 
respect to the potential q{s) but the covariates explain a good part of the 
variability of the potential, then w{s) models only the "residual" random 
field e(s) = q{s) — x(s)/3 and the data y (5) — X {S) (3 can be assumed to be 
not preferentially sampled with respect to e(s). In other words, even when 
the data are preferentially sampled with respect to (?(s), the adoption of 
good (spatial) covariates largely mitigate the problem. If no covariates are 
available and the data are suspected to be preferentially sampled, then the 
approach in Diggle, Menezes and Su (2010) should be considered. 

4.2. Parameter confidence intervals. As known, the classic EM algo- 
rithm does not provide information about the uncertainty of the estimated 
parameter vector ^ . In order to avoid the more cumbersome supplemented 
EM algorithm (see Meng and Rubin (1991)), two methods are proposed to 
solve problem 2. of the above list, namely to provide confidence intervals for 
the elements of ^ . 
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The first method is based on the fact that the maximum likelihood es- 
timator has asymptotically normal distribution A^(^'0)I~^)5 with the 
"true" value of ^ and I the Fisher information matrix. An approximation of 
the information matrix for multivariate normal variables can be evaluated 
as 

(4.9) iij = die'i:-^d,e+^tr{-E-^di-E,^-'dj-E,) 

(see Shumway and Stoffer (2006)), where die and 5jSg are short notations 
for de (*) /d^i and dYl^ (*) /d^i, respectively and 1 <i,j < \^\. 
In the case of the GPM, the vector 

(4.10) e = y-G(l/^ + X/3) 

is normally distributed with variance-covariance matrix 

(4.11) = Far (y - G (1/i + X/3)) 

= Var {--fGw + Ge) 

= G(72Sw + Se)G' 

= gg' (t^Sw + 5:^) 

where is the Hadamard product operator. The solution for the derivatives 
die and SjSg is reported in Appendix B. In the presence of missing data, 
equation (4.9) is still valid but e and Sg have to be replaced with e^^^ and 
Ti^e \ respectively. 

With I available, approximated confidence intervals for the elements of 'I' 
are immediately provided by considering ^^,1"^^ Note, however, that 

N is a good approximation of the distribution \ y (<S)] only when 

N is large, which may not be the case in practical applications. 

To solve this problem, following Fasso and Cameletti (2010), a second 
method based on the bootstrap technique is considered. Let ^ be the es- 
timated parameter vector. For each bootstrap run m, the vector y^^) = 



D 



[m) 



is considered, where 



(4.12) y\'2,=-L-G-(lf, + X0 + ^w^^)+e^ 



and where W(m) and i(m) are realizations from the multivariate normal dis- 
tributions N (O, Se i^e)) ^ ( 0' ( ^) ) 5 respectively. Note that y(^m) 
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preserves the missing data pattern of the observed y. The simulated y(m) is 
used to estimate a new parameter vector ^[m) through the EM algorithm 
and the set 

(4.13) = {l'(i),...,l'(M)} 

is considered as a sample from the distribution | y (5)]. If M is large 
enough, then can be used to derive approximated confidence intervals 
for the elements of ^ without normality assumptions. 

4.3. Potential and conditional potential estimation. Following the plug- 
in approach, the estimated potential is obtained as 

(4.14) g^(s) = /i + X (s) /3 + 7W (s) ; s G P 

where w (s) = (w(s) | y) is the kriging estimate of w (s), which is eval- 
uated analogously to w in equation (4.7). 

The uncertainty of (s) is directly related to the uncertainty of ^ which is 
expressed by | y (5)]. Again, approximated confidence intervals on q^{s) 
can be provided by repeatedly estimating (?5i(s) with ^ extracted either from 

A^^^,I^^^ or from the set defined in equation (4.13). The estimated 
conditional potential is simply given by 

Approximated confidence intervals on q^{s;S) are provided following the 
same approach for q^{s). 

4.4. Total potential estimation. The conditional potential g(s; S) pro- 
vides information about the expected observation when a measuring instru- 
ment is placed at the generic location s given the existence of the other 
instruments. In practice, the following quantity is also of interest: 



(4.15) ^; = max J]]g(s';5\s');5/0 

"^^^ s'e5 

If, for example, q{s) is the spatial market potential, then v represents the 
maximum market volume for the trading area T>. Note that v cannot be 
obtained by simply integrating q{s) or q{s; S) over D. 

A simple way to estimate v is to consider the estimated conditional po- 
tential q^{s]S) with 5 = and to sequentially populate S by choosing the 
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spatial location of T> where g^(s;5) is maximum for the current S. Follow- 
ing this approach, an estimation of v is obtained for \S\ — t- oo. In practice, 
the value of v stops to increase significantly after a finite number of itera- 
tions. Note that a by-product of (4.15) is the optimum S with respect to 
the maximization of v. When the main aim is the optimization of a retail 
network, however, the above approach should be adapted in order to impose 
a threshold on the minimum (geographic) distance between two elements of 
S. 

5. Case study. The GPM is considered here in order to estimate the 
market potential of an economic daily newspaper over the area of the city of 
Bergamo, northern Italy. The aim is to identify the areas with the highest 
market potential, to identify the areas where it would be profitable to open 
additional newsstands and to evaluate the maximum total market volume 
for the city. 

The Italian daily newspaper market is characterized by 64 main newspa- 
per heads with an average market volume of around 5.5 million daily copies. 
As far as the city of Bergamo concerns, only 16 out of 64 newspaper heads 
are commonly commercialized as most of them are local heads referring to 
other Italian cities. The economic newspaper considered in this study rep- 
resents 8% of the total sales volume for the Bergamo area in terms of daily 
copies. Moreover, it should be noted that the economic newspaper is of a 
clientele which differs from that of the most popular newspapers. This im- 
plies that the market potential of the economic newspaper is not necessarily 
reflected in the spatial distribution of the newsstands. In other words, con- 
ditionally to the observed covariates, the sales data are not preferentially 
sampled with respect to the market potential of the newspaper. This is also 
justified by the fact that the sale of daily newspapers represent only 20% of 
the total revenue of a newsstand. 

The data available for the study consist of the yearly average daily number 
of copies sold on working days by = 75 newsstands located over the 
Bergamo area. The sales data of 5 newsstands are unavailable though their 
location is known. The total daily average sales volume for the available 
newsstands is around 491 copies and it is believed that the maximum total 
volume attainable for the city of Bergamo is higher. The newsstand spatial 
locations are shown in Figure 2, along with the circle-plot of the average 
daily number of copies sold. 

By considering the interpretation introduced in Section 3.1, it can be 
stated that the measuring system is represented by the newsstands and 
that the interaction between the measuring instruments is of the absorption 
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Fig 2. Newsstand locations and circle plot of the working day average daily number of 
copies sold. 

type. In fact, once the customer has bought a copy of the newspaper, it 
is absorbed in the sense that the same customer will not buy (during the 
same day) the same copy of the newspaper, neither at the same nor at a 
different newsstand. Since the newspaper price is fixed, the newsstands are 
considered equally-effective and it is supposed that the customer chooses the 
nearest newsstand. Customer loyalty is admitted but it is supposed that a 
newsstand is not more attractive than another. 

The GPM model considered is the same defined in (3.1) but with /i = 0. 
This implies that the market potential goes to zero when moving far from 
the newsstand network as w{s) converges to its marginal mean which is 
zero. It follows that the market potential is zero (or very close to zero) over 
the areas where it would be unfeasible to have a newsstand. The spatial 
correlation function of the latent component w is chosen to be 



while the function (3.2) is considered as the interaction function, with 



(5.1) 




(5.2) 
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Fig 3. Model covariates - (left) spatial density of joint-stock companies; (right) function 
of the geographic distance to the nearest busy street section. 

Table 1 

Estimated model parameter and 95% bootstrap confidence interval 

Pi P2 Sri y e 4> 

estimated 18.29 27.65 11.77 14.63 81.77 231.69 
LCL -0.08 19.19 3.22 11.17 14.92 195.73 

UCL 66.91 66.51 103.33 31.60 253.06 474.64 

Two covariates are considered. The first covariate represents the spa- 
tial density of the joint-stock companies with registered offices in Berg- 
amo. The companies are expected to induce a higher sales volume at the 
near newsstands. The second covariate is a function of the minimum Eu- 
clidean distance to the busiest street sections in terms of people and car 
traffic. In particular, the covariate value at the generic location s is given 
by l/((imm(s) + 0.1), where (imm(s) is the Euclidean distance (expressed in 
kilometers) from the location s to the nearest street section. Both the co- 
variates are depicted in Figure 3. In order to make the /3 coefficients directly 
comparable, the covariates at the newsstand locations have been rescaled to 
the range [0, 1]. 

The model parameter vector ^ is estimated by means of the EM algorithm 
as discussed in Section 4.1 and by using the software provided in [Finazzi 
(2012)]. The estimation result is reported in Table 1 with confidence intervals 
evaluated by following the bootstrap approach discussed in Section 4.2 and 
M = 922. Namely, 95% confidence intervals are obtained by evaluating 
empirical distributions on = •••1 |- Note that the original 

number of bootstrap runs was M = 1000 but 78 runs have been ignored 
after testing the estimated parameters against anomalous values. In fact, the 
EM algorithm is not guaranteed to converge to the global maximum of the 
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likelihood function. In this particular application, the condition (j) > 1500 m 
has been considered to identify bad estimation results as values of (p higher 
than 1500 m implies a very strong and unrealistic competition between the 
newsstands. 

The empirical variance-covariance matrix of ^ is reported in Table 2 and 
it can be compared with the approximated Hessian matrix evaluated by 
considering equation (4.9) and reported in Table 3. In particular, it can be 
noted that the approximated Hessian matrix tends to underestimate the 
variances related to the elements of ^. 

As expected, the (3 coefficients related to the covariates are both positive 
in sign. The coefficient /3i related to the spatial density of the joint-stock 
companies is characterized by a larger confidence interval with lower control 
limit —0.08. Since the confidence interval is an approximation based on 
bootstrap runs, the covariate is retained. 

The estimated 6 ~ 82m suggests that the potential q, net of the covariate, 
is not highly spatially correlated. Moreover, as supported by </> ~ 231m, the 
competition between nearby newsstands is quite strong and two newsstands 
200m apart measure/absorb (on average) only 70% of the actual market 
potential at their locations. 

Given ^, considering equation (4.14), the market potential q^{s) is es- 
timated over the area of the city of Bergamo as depicted in Figure 4. For 
each s £ V, q^{s) provides the daily average newspaper number of copies 
that would be sold by a newsstand if placed at s without any other news- 
stand in v. The maxima of 9|,(s) correspond to commercially strategic areas 
that should be served by at least one newsstand. The global maximum is 
equal to 79.86 yearly average newspaper copies and it is located at 45.6930° 
latitude and 9.6640° longitude. The estimated latent variable w{s) is dis- 
played in Figure 5. It can be noted that its role is more pronounced in the 
city center where, apparently, the covariates are less capable of explaining 
the observed market potential. The estimated conditional market potential 
q^{s; S), depicted in Figure 6, provides the daily average newspaper number 

Table 2 

Empirical variance-covariance matrix of 'if based on 922 bootstrap runs 





/3i 


/32 


ol 


7 


e 


</> 


Pi 


303.00 


31.87 


207.95 


50.04 


108.58 


700.60 






217.28 


450.01 


90.42 


899.77 


1091.76 


ol 






6180.13 


206.79 


4343.86 


2725.23 


7 
B 










51.04 


30.94 
8941.41 


578.37 
693.87 
7268.92 



GEOSTATISTICAL MODELING UNDER SPATIAL INTERACTION 17 




Longitude 

Fig 4. Estimated potential q^{s) (average daily number of copies) over the area of the city 
of Bergamo. 

of copies that would be sold by a newsstand if placed at s given the current 
newsstands located at S. Thus, the maxima of q^{s; S) represent the spatial 
locations where it would be profitable to open additional newsstands. Fi- 
nally, the bootstrap standard deviation of the estimated conditional market 
potential, representing its uncertainty, is shown in Figure 7. 

The total market volume related to the economic newspaper and the Berg- 
amo area has been evaluated following the procedure discussed in paragraph 
4.4. Figure 8 shows the total market volume and its 95% confidence interval 
with respect to the number of newsstands. Note that the market volume 

Table 3 

Approximated Hessian matrix for ^ 
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Fig 5. Estimated latent variable ■w{s) over the area of the city of Bergamo. 



stabilizes at around 688 average daily copies after 200 newsstands. This is a 
consequence of the fact that 200 newsstands absorb most of the market po- 
tential and add more newsstands does not increase the total market volume 
significantly. Also note that the optimized retail network of 75 newsstands 
absorbs a market volume equal to 646 copies, which corresponds to 94% of 
the maximum market volume and is 30.8% higher than the market volume 
absorbed by the actual retail network. 

In the light of this results, the publisher of the economic daily newspaper 
can consider to improve the retail network in order to increase the daily 
market volume. Additional newsstands can be opened only with the consent 
of the municipal authority and, since the economic newspaper represents a 
small part of the daily revenue of a newsstand, it cannot be guaranteed that 
the additional newsstands would be opened where the market potential of 
the newspaper is high. Nevertheless, the Italian market of daily newspapers 
is undergoing a liberalization phase and the publisher should start thinking 
new forms of retailing. 

6. Conclusions. The geostatistical potential model has been proven to 
be an essential statistical tool for the estimation of the spatial market poten- 
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Fig 6. Estimated conditional potential q^{s;S) (average daily number of copies) over the 
area of the city of Bergamo. 

tial of a retail product from its sales data. The model output is immediately 
and easily interpretable (uncertainty included) as it is provided in the form 
of spatially continuous surfaces and with the same unit of measure of the 
original data. 

The geostatistical potential model has been successfully applied to the 
estimation of the spatial market potential and to the total market volume 
of an economic daily newspaper for the city of Bergamo, Italy. The analysis 
of the results allowed to conclude that the daily sales volume can be signifi- 
cantly increased by focusing on the areas of the city which are characterized 
by a high market potential but they are not properly covered by the retail 
network. 

Future extensions of the model include the introduction of the time vari- 
able, in order to describe and study the temporal fluctuations of the spatial 
market potential, and the relaxation of the equally-effectiveness property, 
in order to address the case of stores characterized by a different degree of 
attractiveness. 

As a final remark, it is worth noting that the geostatistical potential model 
can be applied outside the geomarketing field. For instance, the sales data 
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Longitude 

Fig 7. Bootstrap standard deviation map of the estimated conditional potential q^{s;S). 

of the chemists of a city with respect to a drug can be analyzed in order to 
assess the diffusion of a disease in terms of a spatially continuous surface. 
More generally, the model can be applied in the case where the data related 
to a set of statistical units are available in aggregated form (e.g. due to 
privacy reasons) but they are georeferenced with respect to precise points 
in space. 

APPENDIX A: LATENT VARIABLE ESTIMATION 

The Gaussian latent variable w is estimated by applying the usual for- 
mulas of the multivariate normal distribution. In particular: 

w = E^(k) (w I y) 

(A.l) = Swyl^yMy - G (1^ + X^)] 

where 

Sy = Far [G + X/3 + 7w + £)] 

= GVar[-fw + e]G' 

= G (t^Ew + Ss) G' 
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Fig 8. Total market volume with respect to the number of newsstands and 95% confidence 
interval based on 922 bootstrap runs. 



and 

Swy = E [{w - 0) ■ [y - E {y)]'] 

= E [w (7Gw)'] 

= l^^G' 

The variance of the estimated w is given by 

A = Var^^k) (w | y) 



(A.2) 



^wySy ^ (^wy)' 



When y is characterized by missing data, equations (A.1-A.2) become 



w 
A 



[SwyL') (LSyL') ML(y-G(l^ + X/3))] 



wyy 



APPENDIX B: VECTOR AND MATRIX DERIVATIVES 

The evakiation of the approximate Fisher information matrix defined in 
equation (4.9) requires the computation of the vector derivatives de (^') /d^i 
and the matrix derivatives SSg (^) /d^i, 1 < i < l^*], with e and 'S^ defined 
in equations (4.10) and 4.11), respectively. 
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In the case of the spatial correlation function defined in equation (5.1) and 
the interaction function defined in equation (5.2), the following derivatives 
hold: 



(B.l) 



(B.2) 



d-^e (^) 



-g if *i = /u 

-gQx; if '^i = fSi; l<l <b 

-ag^ (1/x + X/3) if ^, = ^ 

otherwise 

gg' In if *i = (Te 

2jgg' if = 7 

T'gg'oioSw if '^^ = 



G (t^Sw + S,) 




if ^z = (p 
otherwise 



where is the / — th column of the matrix X and H is the distance matrix 
based on S. 

Finally, the {pq) — th element of the matrix G is given by dgp ■ gq + Qpdgq , 
where gp is the p — th element of the vector g and while the p — th element 
of the vector (9g^ is given by 



(B.3) 



dgp 



dgp 



exp 



N 

^ exp 



with hpq the {pq) — th element of the matrix H. 

SUPPLEMENTARY MATERIAL 

Supplement A: Dataset and Matlab@ code 

(http://lib.stat.cmu.edu/aoas/). Georeferentiated newsstand sales data and 
Matlab@ code for the data analysis. 
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